#pragma once

#include <glib/g2/smatrix.h>
#include "outline/type_traits.h"

namespace outline {

/// fixed_line_estimator provides a caller with discrepancy of the line, that must pass through the first and the last points,
/// from the point set; incremental modification of the point set is allowed.
///
/// The class accumulates information about start, several intermediate, and finish points,
/// and provides summed square of discrepancy of intermediate points from line start-finish.
/// The class allows easy testing of new finish point (previous finish point becomes else one intermediate point).
template <typename T>
class fixed_line_estimator;

template <typename T>
class fixed_line_estimator< point<T> > {
public:
	typedef point<T> sample;
	struct patch {};
public:
	typedef typename type_traits<T>::accumulator_type accumulator_type;
	typedef g2::basic_vector<accumulator_type>	vector;
	typedef g2::basic_smatrix<accumulator_type>	smatrix;
public:
	fixed_line_estimator() { }
	fixed_line_estimator(vector const & s)
		{ start(s); }
	static int min_step_len()
		{ return 2; }
public: // methods, required by simplifier
	///forgets accumulated information and sets new start point of line
	void start(vector const & s)
		{ r0 = s; q.zero(); }
	///accumulates information about else one intermediate point
	void add_sample(vector const & i)
		{ q += tsquare(i - r0); }
	///accumulates one point and returns discrepancy,
	///optional argument may return parameters of the line (reserved for objects other from lines)
	double add_sample_and_estimate(vector const & i, patch * = 0, bool swap_ends = false) {
		smatrix sqr_delta = tsquare(i - r0);
		double r = measure(sqr_delta);
		//assert (r >= min_discrepancy()); //may be not satisfied due to rounding errors
		q += sqr_delta;
		return r;
	}
	///same as add_sample_and_estimate(...) but adding point coincides with start point
	///typically this method is called the latest for closed contours
	double add_start_and_estimate(vector const & i, patch * = 0, bool swap_ends = false) {
		assert (i == r0);
		return discrepancy_from_start();
	}
	///wherever possible added and final point are, discrepancy will not be less than returned value
	double min_discrepancy() const
		{ return g2::smatrix(q).eigenvalue(0); }
public: // additional methods, not called from simplifier
	///what discrepancy would be if finish is in given point?
	double finish(vector const & f) const
		{ return measure(tsquare(f - r0)); }
	///measures distances from start point (not from any line)
	double discrepancy_from_start() const
		{ return q.tr(); }
private:
	double measure(smatrix const & sqr_delta) const
		{ return q.tr() - double(dscalar_product(q, sqr_delta)) / sqr_delta.tr(); }
private:
	vector r0;
	smatrix q;
};

/// best lines for all the channels
template <typename T>
class fixed_line_estimator< std::vector< point<T> > > {
public:
	typedef std::vector< point<T> > sample;
	struct patch {};
public:
	fixed_line_estimator() { }
	fixed_line_estimator(sample const & s)
		{ start(s); }
	static int min_step_len()
		{ return 2; }
public:
	void start(sample const & s) {
		channels.resize(s.size());
		for (size_t i = 0; i < channels.size(); ++i)
			channels[i].start(s[i]);
	}
	void add_sample(sample const & s) {
		assert (s.size() == channels.size());
		for (size_t i = 0; i < channels.size(); ++i)
			channels[i].add_sample(s[i]);
	}
	double add_sample_and_estimate(sample const & f, patch * = 0, bool swap_ends = false) {
		assert (f.size() == channels.size());
		double m = 0;
		for (size_t i = 0; i < channels.size(); ++i)
			m += channels[i].add_sample_and_estimate(f[i]);
		return m;
	}
	double add_start_and_estimate(sample const & f, patch * = 0, bool swap_ends = false) {
		assert (f.size() == channels.size());
		double m = 0;
		for (size_t i = 0; i < channels.size(); ++i)
			m += channels[i].add_start_and_estimate(f[i]);
		return m;
	}
	double min_discrepancy() const {
		double m = 0;
		for (size_t i = 0; i < channels.size(); ++i)
			m += channels[i].min_discrepancy();
		return m;
	}
public:
	double finish(sample const & f) const {
		assert (f.size() == channels.size());
		double m = 0;
		for (size_t i = 0; i < channels.size(); ++i)
			m += channels[i].finish(f[i]);
		return m;
	}
	double discrepancy_from_start() const {
		double m = 0;
		for (size_t i = 0; i < channels.size(); ++i)
			m += channels[i].discrepancy_from_start();
		return m;
	}
private:
	std::vector< fixed_line_estimator< point<T> > > channels;
};

} //namespace outline
